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We consider the distribution of the transmission coefficients, i.e. the singular 
values of the modal transmission matrix, for 2D random media with periodic 
boundary conditions composed of a large number of point-like nonabsorbing 
scatterers. The scatterers are placed at random locations in the medium and 
have random refractive indices that are drawn from an arbitrary, known 
distribution. We construct a randomized model for the scattering matrix that 
retains scatterer dependent properties essential to reproduce the transmission 
coefficient distribution and analytically characterize the distribution of this 
matrix as a function of the refractive index distribution, the number of 
modes, and the number of scatterers. We show that the derived distribution 
agrees remarkably well with results obtained using a numerically rigorous 
spectrally accurate simulation. Analysis of the derived distribution provides 
the strongest principled justihcation yet of why we should expect perfect 
transmission in such random media regardless of the refractive index distribu¬ 
tion of the constituent scatterers. The analysis suggests a sparsity condition 
under which random media will exhibit a perfect transmission-supporting 
universal transmission coefficient distribution in the deep medium limit. @ 
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Introduction 


Materials such as turbid water, white paint, and egg shells are considered opaque because 
multiple scattering by the randomly placed constituent scatterers in the medium frustrates 
the passage of light . The seminal papers by Dorokhov , Barnes and Pendry et al. [3,18 


and others [4,17 postulate that even if a normally incident wavefront barely propagates 
through a thick slab of such media, there will generically exist a few highly-transmitting 
wavefronts that will propagate through the slab with a transmission coefficient close to 1, 
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i. e, they will be nearly perfectly-transmitting. These perfectly transmitting eigen-wavefronts 
are the right singular vectors of the modal transmission matrix and are optimized to the 
specihc random medium. 

These seminal papers inspired the breakthrough experiments by Vellekoop and Mosk 
27 , and others iillzll 14, [^[^[^1^1^ provide credence to the hypothesis that there 


generally exist (nearly) perfectly transmitting eigen-wavefronts in highly scattering random 
media composed of a larger number of non-absorbing scatterers. Recently, we verihed this 
hypothesis 12,13 for 2-D systems with periodic boundary conditions composed of hundreds 
of thousands of non-absorbing scattering using numerically rigorous simulations. 

The perfect-transmission supporting universal transmission coefficient distribution postu¬ 
lated by Dorokhov, Mello, Pereyra, Kumar, Pendry, and Barnes 


17,18 , was derived 


assuming that the medium was deep enough so that the scattering matrix obeyed a phys¬ 
ically consistent {i.e. obeying reciprocity and time-reversal conditions) maximum-entropy 
law. Their analysis does not provide a principled and mathematically grounded framework 
for reasoning about whether, when or the sense in which a deep medium composed of a large 
number of randomly placed point-like scatterers with an arbitrary distribution of refractive 
indices can be expected to have a perfect transmission-supporting transmission coefficient 
distribution. 

In this paper, we use modern random matrix theory to revisit the problem of predicting 
the transmission coefficient distribution of 2-D random media with periodic boundary condi¬ 
tions composed of a larger number of randomly placed point-like scatterers with an arbitrary 
refractive index distribution. We provide a characterization of the transmission coefficient 
distribution that explicitly depends on the refractive index distribution, the number of prop¬ 
agating modes and the depth of the medium for layered random media (in a sense we will 
make precise) composed of a large number of point-like scatterers. 

The critical part of our derivation relies on the development of an isotropic random matrix 
model for the modal transfer matrix of a single randomly placed point-like scatterer. The 
random transfer matrix has singular value distribution that matches the singular values of 
the physical transfer matrix of a randomly-placed point-like scatterer. However, the left and 
right singular vectors of our random transfer matrix construction are modeled as indepen¬ 
dent and isotropically random. This allows us to use tools from free probability theory to 
approximate the transmission coefficient distribution of a layered random media composed 
of layers containing point-like scatterers. 

We show that the derived distribution agrees remarkably well with results obtained us¬ 
ing a numerically rigorous spectrally convergent simulation that utilizes spectrally accurate 
methodologies. This justihes the use of our isotropic model for reasoning about the properties 
of the derived distribution. Analysis of the resulting distribution brings into sharp focus the 
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universal, i.e., scatterer-property independent, aspects of the distribution and provides the 
strongest principled justihcation yet of why we should expect perfect transmission in such 
deep random media regardless of the refractive index distribution of the constituent scat- 
terers. The analysis brings into focus a sparsity condition under which random media can 
be expected to exhibit a perfect transmission-supporting universal transmission coefficient 
distribution in the deep medium limit. 

We describe the setup and dehne the transmission coefficient distribution in Section]^. We 
highlight some pertinent properties of the system modal transfer matrix in Section and 
employ them in Section]^ to formulate a isotropically random model for the transfer matrix 
of a single point-like scatterer. In Section we describe the pertinent free probabilistic 
tools from random matrix theory that allow us to analytically characterize the limiting 
transmission coefficient distribution of a medium composed of many scatterers from the 
eigen-distribution of the isotropic transfer matrix of a single point-like scatterer. We analyze 
the properties of the limiting transmission coefficient distribution thus obtained in Section 
and bring into sharp focus its universal, i.e., scatterer property independent, aspects. 
We validate our theoretical predictions using numerically rigorous simulations in Section 
Details of some computations have been relegated to the Appendix. 

2. Setup 

We study scattering from a two-dimensional (2D) random slab of thickness L and periodicity 
D; the slab’s unit cell occupies the space 0 < x < D and 0 < y < L (Fig.[^. The slab contains 
A’lay inhnite and z-invariant circular cylinders of radius r that are placed randomly within 
the cell, as described shortly. The cylinders are assumed to be dielectric with refractive index 
Ud', care is taken to ensure the cylinders do not overlap. The radius of the cylinders is chosen 
to be much smaller than the wavelength A so they, in effect, act like point scatterers. 

For ic = 1,2,..., N\ay, the x and y position of the center of the ic-th cylinder are 
Ux,ia,Uy,ic + (A — 1)-^), respectively where and Uy^i^ s are i.i.d, random variables with 
uniform distribution on [r,D — r] and [r,i — r], respectively. Here £ = L/N\ay is depth of each 
“layer”; £ is chosen to be larger than y/DX. Each cylinder’s refractive index Ui^ is drawn 
independently from the same distribution of refractive indices ri{n). 

Fields are TM^ polarized: electric helds in the ?/ < 0 (i = 1) and y > L = N\ 3 y£ {i = 2) 
halfspaces are denoted ej(p) = ei{p)z. The held amplitude ei{p) can be decomposed in terms 
of +y and —y propagating waves as ei{p) = e'^{p) + e^{p), where 

N 

4(p)= (1) 

n=-N 

In the above expression, p = xx + yy = {x,y), = kn,xX ± kn,yy = {kn,x,±kn,y), kn,x = 


3 


Periodic repitition 
with period D 



D 


“2 


“2 


Fig. 1. Setup. 


2Tm/D, kn,y = 27r^(l/A)2 - {n/DY, A is the wavelength, and = \J\\]^\\ 2 /kn,y is a 
power-normalizing coefficient; a time dependence is assumed and suppressed. We assume 
N = [D/AJ, i.e. we only model propagating waves and denote M = 2N + 1. The modal 
coefficients i = 1,2] n = —N,..., N are related by the scattering matrix 
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where = af_p^ ... af^ ... af^ . In what follows, we assume that the slab is only 
excited from the y < 0 halfspace; hence, = 0. For a given incident held amplitude ei{p), 
we dehne the transmission coefficient as 


T[a{) := 


\\S2i-at\\l 


ll«^ 


+ 112 
2 


(3) 


We denote the transmission coefficient of a normally incident wavefront by Tnormai = 

1 T 


); here ^ denotes transposition. 
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2.A. The transmission coefficient distribution 

The problem of designing an incident wavefront that maximizes the transmitted power 
can be stated as 

/ +\ 11*^21 ' tki II 2 no +ii2 / A\ 

= argmax r= arg max — ^ — = arg max ||^21 ■ at ||2 

a+ a+ il% II2 l|ai''l|2=l 

where || atf || 2 = 1 represents an incident power constraint. 

Let S '21 = ■ V 4 denote the singular value decomposition (SVD) of S' 21 ; cxj is the 

singular value associated with the left and right singular vectors and v^, respectively. By 
convention, the singular values are arranged so that ai > ... > au and ^ denotes complex 
conjugate transpose. Then via a well-known result for the variational characterization of the 
largest right singular vector Theorem 7.3.10] we have that 


“opt = hi- 


(5) 


When the optimal wavefront a^p^ is excited, the transmitted power is Xopt := = af. 

When the wavefront associated with the Tth right singular vector is transmitted, the 
transmitted power is r* := T{vf) = af, which we refer to as the transmission coefficient of 
the Tth eigen-wavefront of S' 21 . We are interested in the limiting transmission coefficient 
distribution whose p.d.f. is dehned as 


fir) = lim E 

M,Nt^y^oa 


M 
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^h(r-r(nj) 


2 = 1 


lim E 

A^iA^lay— ^CX) 


M 


M 






2=1 


( 6 ) 


where we assume that A^iay/M —)■ c G (0, cxd) as M, iV|ay —)■ 00. The Dorokhov-Mello-Pereyra- 
Kumar (henceforth, DMPK) distribution [8,17 has density given by 


/DMPK('r) — 


free 


2L T\/l — r ’ 


for 4exp(-L/2/ 


free. 


< 


'T < 1 , 


(7) 


where /free is the mean-free path in the medium. The DMPK distribution is posited 


17,18 to be the universal limiting distribution for systems comprised of many scatterers in 


the limit where L ^ M. 

Assuming a scattering regime where the DMPK distribution holds, Eq. ([^ predicts the 
existence of highly-transmitting eigen-wavefronts that achieve (nearly) perfect transmission. 
Since the DMPK distribution was derived under a maximum-entropy type assumption (which 
we shall revisit shortly), the material properties of the scatterers, such as the distribution 
of refractive indices, do not explicitly appear in the expression in Eq. ([^ for its p.d.f. but 
instead are encoded implicitly via the /free parameter. Our objective is to theoretically predict 
/(r) in Eq. ([^ and explicitly characterize its dependence on the refractive index distribution 
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r]{n), N\s,y, and M, assuming we are in a regime where each scatterer is small enough so that 
it effectively acts as an isotropic point scatterer. Our mathematically-derived framework 
permits reasoning about the conditions under which we might expect a universal limiting 
distribution and the existence of the (nearly) perfectly transmitting eigen-wavefronts. 


3. Background: the transfer matrix and its pertinent properties 

The scattering matrix S in Eq. ([^ describes the relationship between the modal coefficients of 
incoming and outgoing waves. Rearranging the terms in Eq. (|^ relates the modal coefficients 
in i = 1 and i = 2 halfspace via the transfer matrix T 


at 


c c c 

*^21 “ *^22 * ^12 ' *^11 

1 

d 

(M 


at 

.—2 . 

'S. 

-5i 2-' ■ 5n 

0-1 

^^12 J 

✓ 

-1 

cil. 


=:T 
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where we have assumed that the 512 matrix is invertible. Rewriting the transfer matrix as 


T = 


S22 

1 

0 

- 1 

1- 

0 

I 

d 



Q —1 Q Q~^ T 

^22 * *^21 ‘ *^11 *^12 ^ 

-I 


1 - 

1 

1 - 

0 

- 1 

1 

0 

0-1 

^^12 


( 9 ) 


allows us to easily verify that det(T) = det(T^-T) = 1. In the lossless setting when -S = I, 
and 5 i 2 is invertible, it is shown in Appendix that the 2M eigenvalues of ■ T denoted 
by Ai > ... > \ 2 m are 

2-Ti + 2^Jl - n 2-Ti- 2V1 - Ti t ■ ^ 

Ai = - and A2M-i+i = -y- tor ? = 1,..., M. (lU) 


Ti 


Ti 


Note that Ai ■ A 2 M-i+i = 1 so that the 2M eigenvalues of ■ T come in reciprocal pairs. 
From Eq. (10), we have that 


Ai + \ 2 M-i+l — -2, 

Ti 


so that 


Ti = 


Ai + 1/Ai -|- 2 
Substituting A* = exp(2xi) in Eq. (0 yields 

4 1 


( 11 ) 


Ti = 


exp(2a;i) exp(-2a:i) -h 2 (exp(xi) exp(-a:i)/2)2 cosh^(a;i) 
Equivalently, since Xi = 0.5 InAi, we have 

1 


Ti = 


cosh^(0.5 In Ai) 


fA Ai = exp(2cosh AI/^/t^)), 


( 12 ) 
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and we have obtained a direct relationship between the eigenvalnes of -T and the transmis¬ 
sion coefficients. Let h{X) denote the limiting eigenvalne distribntion of the transfer matrix 
dehned as 

.. 2 M 


h(A) = lim E 

A</,7V|3y—^OO 


2M 


i=l 


(13) 


Then, a direct conseqnence of Eq. (12) is that once we know h(A), a simple change of variables 
yields the transmission coefficient distribntion /(r) as 


/(r) = h(A) 


\dT/dX\ 


A in Eq. | [T^ 


A=exp(2cosh ^{l/y/r)) 


(14) 


Since the eigenvalnes of ■ T come in reciprocal pairs, h(A) for A < 1 nniqnely determines 
h{X) for A > 1. Thns we can rewrite Eq. (14) as 

(A+ 1)3 


/(r) =2h(A) Ia<i 


4|A 


A=exp(2cosh ^(l/y/r)) 


(15) 


where Ia<i denotes the indicator fnnction on the set A < 1. From Eq. (0 , we have that 

(A + 1)^ 1 _(A + 1)^ 


1 

r 


4A 


and 


r 


(A-1)2’ 


so that rearranging terms on the right hand side of Eq. (|15|), yields 

1 


/M = 


Ta/I — T 


■{2h(A)AlA<i} 


(16) 


A=exp(2cosh ^(1/^^)) 


We note that Eq. (16) is an exact relationship between the eigenvalue distribution of the 
transfer matrix and the transmission coefficient distribntion. Comparing Eqs. (16) and ([^ 
reveals the important insight that the DMPK distribntion arises nnder the assnmption that in 
the limit of deep random media, h{X) = /free/(4LA), or eqnivalently, that h{X) is a log-uniform 
distribution. This is the maximum-entropy assumption that yields the DMPK distribntion 
for deep random media. Onr goal is to analytically characterize h{X) and hence /(r), via Eq. 
(16) as a fnnction of iV|ay, M and the refractive index distribntion of the scatterers for the 
setnp in Fig. 


4. An isotropically random model for the transfer matrix of a single point-like 
scatterer 

Let Tj denote the transfer matrix of a layer containing a single scatterer (Fig. and let 
*S'u , S 22 , and S 21 denote its scattering matrix and snbblocks thereof, respectively. When the 
scatterers are point-like and D is large, then S*!! and S '22 are well approximated by a rank 
one matrix whose largest singnlar valne a G [0,1) we will refer to as the scattering strength. 
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This is obviously true for D —)■ oo, and remains remarkably accurate for smaller D as well. 
Since S is unitary for lossless media, we have that ■ Sn + S 21 ■ S 21 = I- Hence, the 521 
matrix must have an SVD of the form 


5® = U, ■ diag(l,..., 1, 

where f/j and Vi are the left and right singular vectors of which encode the physics 
of the scattering system. Consequently, by Eq. §, the transmission coefficients of S 2 I are 
approximately 

Ti ~ ... ~ tm -1 ~ 1, and tm ~ 1 - 

From Eq. (10), we can conclude that the 2M — 2 eigenvalues of ■ T will equal one. The 
remaining two eigenvalues will Am and Am+i = 1/ Am where 


_ 2-(1-q^) + 2^1-(1-q2) _ 1 + + 2 g _ 1 + g 

^ 1 — g^ (1 —g)(l + g) 1 —g 

This implies that the transfer matrix will have an SVD of the form 

T, = U, diag(l, 1, • • ■ , 1,1, l/Ve) 

"V" ^ 

2M-2 entries 


(17) 


( 18 ) 


where Ui and V are the left and right singular vectors of Ti, which again encode the physics 
of the scattering systems. The refractive index distribution 7 ]{n) induces a distribution feit) 
on 6 which we assume to known and obtained either a via a change of variables as 


g 


V 16D 



(19) 


under a point scatterer assumption for the large D, r A, and rid ^ 1 regime, or using 
computational electromagnetic techniques. 

The transfer matrix of the entire system in Fig. is obtained from those of the layers as 


Af|ay 

t = 1[t, 

i=l 


( 20 ) 


Each of the transfer matrices are independent and identically distributed. Fig. plots the 
expected values of the squared magnitude of the (bistochastic) correlation matrix formed 
by the inner product of the right singular vectors of a transfer matrix associated with a 
single randomly placed scatterer and the left singular vectors of an independent transfer 
matrix associated with another randomly placed scatter, averaged over 100,000 independent 
realizations. If the singular vectors were independent and isotropically random (or Haar 
distributed) then we would get an empirically averaged matrix with all of its entries close 













to 1/M. From Fig. we can conclude that the singular vectors of two independent transfer 
matrices are not isotropically random with respect to each other. However, most of the entries 
of the correlation matrix have entries ‘close’ to 1/M. 

Very recently, Anderson and Farrell rigorously showed that the product of independent 
(Hermitian) random matrices with independent eigenvectors having a correlation matrix 
whose entries have squared magnitude entries exactly equal to 1/M will have the same 
limiting distribution as the product of independent random matrices with the same eigenvalue 
distribution but isotropically random eigenvectors. Here too, we have a situation where we 
are interested in analyzing the singular value distribution of products of random matrices 
with independent singular vectors. However, the correlation matrix of the singular vectors 
has entries whose squared magnitude is not exactly 1/M, as would be the case if the left 
and right singular vectors were isotropically random, but instead close to 1/M. This leads 
to our conjecture that the singular value distribution of the matrix in Eq. (20) can be ‘well 
approximated by’ the singular value distribution of independent random matrices with the 
same per-matrix singular value distribution but isotropically random left and right singular 
vectors. 

Motivated by this conjecture, we now consider an isotropic random matrix model for the 
transfer matrices whose singular values are specified by Eq. (18) but whose left and right 
singular vectors are independent and isotropically random. We then use tools from free 
probability theory to analytically characterize the transmission coefficient distribution that 
arises due to this isotropic model for the transfer matrix of a point-like scatterer. Numerically 
rigorous physical simulations in Section will validate our conjecture. A mathematically 
rigorous treatment of this conjecture, including a quantihcation of the approximation error, 
remains an open problem. 
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(a) 3D plot. (b) Top view. 


Fig. 2. Relationship between singnlar vectors plotted in lOlog^Q scale. The 
absolnte valne sqnared of the correlation matrix between the singular vectors 
averaged over 100000 trials. The settings were n = 3.1, a = 0.9, 6 = 0.053, r = 
0.05A, I = 11.66A, D = 50.43A, M = 101. 


5. Analytically characterizing the transmission coefficient distribution 


We now discuss some preliminaries required to compute the transmission coefficient distribu¬ 
tion under the isotropic transfer matrix assumption. Let Xm be an M x M symmetric (or Her- 
mitian) random matrix whose ordered eigenvalues are denoted by Ai(Xm) > • • ■ > Xu^Xm)- 
Let hxj^f be the empirical eigenvalue distribution, i.e. , the probability distribution with 
density 

1 ^ 
i=l 

Now suppose that Am and Bm are two independent M x M matrices whose empirical 
eigenvalue distributions converge as M —> oo to non-random distributions having densities 
Ha and respectively. A natural question then is: how is the limiting eigenvalue distribution 
of the matrix ■ Am ■ Bm related to the limiting eigenvalue distributions of Am and 


Sm? 


Free probability theory 28,29 states that if we know Ha and hs and the matrices Am and 
Bm are asymptotically free, we can compute the limiting eigenvalue distribution of Am ■ Bm 
from the limiting eigenvalue distributions of A and B. Specifically, in this setting, Iia-b is 
given by the free multiplicative convolution of Iia and Iib, denoted by Iia^ Iib which is 
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computed as described next. We first define the S -transforn^ which is given by 

1 + ^ 1 


^I’xiz) := 




where 


and 


^x{z) - 


z — t 


9x{z) = 


hx{t)dt = -1 + zgx{z), 


hxit)dt, 


z — t 


is the Cauchy transform of hx- Then S'-transform of hx-B is 

'4’A-Biz) = ^pAiz)^jJB{z). 


( 21 ) 

( 22 ) 

(23) 

(24) 


Note, that given the Cauchy transform gx{z), we can recover the density via the inversion 
formula ^ 

hx (z) = -lim Im gx (z + j e). (25) 

TT e-S>0 

A sufficient condition for the asymptotic freeness of two random matrices is that their sin- 
gnlar vectors are independent and isotropically random . Consequently, under the isotropic 
transfer matrix assnmption, the transfer matrices of snccessive layers are asymptotically free, 
by construction. Hence, we can use free multiplicative convolution machinery to characterize 
the limiting eigenvalne distribution of the transfer matrix of a mnlti-layered scattering system 
as depicted in Fig. since, by Eq. (20), the composite transfer matrix is the product of A"iay 
independent (and asymptotically free) random transfer matrices each having independent, 
isotropically random left and right singnlar vectors and singnlar valnes given by Eq. (18). 

To that end, we hrst compnte the empirical eigenvalue distribution of TT . Tj which is 


h*(A) = 1 


_ , 5(A - 1) + -d) + - 1/0). 

2M ^ ’ 2M ^ ’ 2M ^ ’ 


Its Cauchy transform is given by 


and 


gAz) = 1 — — 

’ ' M/ z 


1 1 
— + 


2M 


z-i/e 


(26) 


(27) 


ii{.z) — —1 + ( 1 — 

= -1 + 


Z 1 

+ 




+ 




M J z — 1 ' 2M \z — 0 z — 1/6 
z 1 ( z 0.5 z 0.5 z 


z-1 
1 1 


M \z-l 


0.5z 0.5z 

+ 


z-6 z-1/6 

z 


( 28 ) 


M\z-e z-i/e z-i 

=AoC) 


■ V' 

=T(d 


^Denoted here by 'i/’(’) to avoid any confusion with the S (or scattering) matrix. 
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Repeating the computation for the setting where the ^j’s are random with pdf fg{-) yields 




0.5 0.5 

+ 


fe{t)dt. 


(29) 


z — t z — 1/t z — 1 

To compute f/’ilz) using Eq. (21) we need to compute A standard application of 

perturbation theory (see, e.g., 1^) yields 


(^) = ^0 (^) - T7 


1 ^(a^) 


M da:^o{x) 


X=^g (Z) 




Substituting ^(z) = [z + 1)/z gives 


= 


z + 1 1 i{x) 


z M d^^o{x) 


+ 0 ( 


x=l+l/z 


M2 J 


or equivalently 


^ + 1 


so that by Eq. (21), 


^ = 1 + — ^ + 0 (—^ 

^ ^ ^ M z + 1 z^ ^ \M^ J 


Mz) = 1 - + O 


Mz{z + 1) \My 


Then, 


A'lay 

i’hiz) = Yli’iiz) = 


2=1 


U(i + i) 

Mz{z + 1) ' \M^ 


+ 0 — 


1 Af|a. 


In the regime where M, A'^iay —?• oo with iV|ay/M —)■ c we obtain 


A{z]c) = lim 

M^oo 


exp 


1 _ +0 f—^ 

Mz{z + 1) \My 


M 


e(i + ^) ' 

z{z + l) 


= exp —c • 


^(1 + 1) ' 

2(2 +1) 


(30) 


(31) 


(32) 


(33) 


We next discuss how to obtain the distribution from 'iph{z;c). Inserting z = ^h{y) in Eq- 

(34) 


(21), we get 


1 + ^h{y) 1 / /'c z' 

-- = ^h{^h{y))- 


^h{y) y 


Substituting in the expression for ihiy) from Eq. (22), we get 

1 - ^ + ygh{y) 1 , / . , ^33 

—y—-yx— = 'il^hi.-l + yghi.y))- 

-^ + ygh[y) y 


(35) 
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Therefore, we get the fixed-point equation 


9h{z) 


zgh{z) 


= 'tphizghiz) - 1 ). 


Substituting Eq. (33) 


9h{z) 


Z9h{z) - 1 


= exp 


e( 


—c ■ 


zghjz) ' 

zgh{z)-i> 


(:3b(~) + 1) 


(36) 


The density h{X) can be recovered from the Cauchy transform gh{X) using Eq. (25) after 


solving the hxed-point equation. The transmission coefficient distribution is then obtained by 


Eq. (16). Note that .^(^) in Eq. (29) explicitly encodes the portion of the limiting distribution 


that depends on the scatterer-dependent properties via fe(t), where 9 is related to the scat¬ 


tering strength a of a single scatterer via Eq. (17) and a is related to the scatterer-dependent 


properties via Eq. (19). 


6. Properties of the limiting transmission coefficient distribution 


We now analyze the properties of the distributions characterized by Eq. (36). The mean of 
f{r) is 

E[r] = [ Tf{T)dT = [ , , , h{X)dX = 4 [ ^. h{X)dX (37) 


A + A-i + 2 


(A + 1)^ 


where we have used Eq. (11) to express E[r] with respect to h{X). From Eq. (22), we note 
that 


:= d,^h{z) = - 


X 


-h{X)dX. 


Thus by comparing the righthand sides of Eqs. (37) and (p^, we have that 


E[t] = -45i(-l). 

From the computation in Appendix we obtain the closed-form expression 

E|rl =- ^ ^ 


(38) 


(39) 


l + c 


t 


1 -l- t 


fe{t)dt 


1 + cBo 


(40) 


= .B2 


Here, we call B 2 the normalizing factor, and it represents the average scattering strength of 
a single layer. The normalizing factor can be used to homogenize two different materials by 
giving measures to calculate the effective lengths, and its specific usage will be discussed in 
Section]^ We now compute the second moment of /(r), which is given by 

16 ... /■ A^ 


E[r ] = / r f{T)dT = 


(A + A-i + 2)2 


h{X)dX = 16 


(A + 1) 


-h{X)dX. (41) 
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Note that 


and 

so that 
1 


ga"(^) - = 


C{z) := d,i'j^{z) = -2 

C(^) := d^a{z) = -6 


A A 

+ 


A 


(A-;.)3 

A 

(A^ 


h(A)dA, 

■h(A)dA, 


(A — zY (A — zY 


h{X)dX = 


X^-Xz-X 


h{X)dX. (42) 


Comparing Eqs. (43) and (42) gives ns the relationship 

'1 


EM = 16 [^C(-l) - ^a'(-l)J ■ («) 

The closed-from expression for the second moment is lengthy and derived in Appendix [Cj 


From Eq. (67) and Eq. (60) we obtain 


E[r=] 16 - §«(-!) 

E[r] -4 


a(-l) 


(44a) 


-iV-i'" 

h S/x 


6 


= -4- 




/X 5 


+ 1 


- 1 " 


2 (C") 


/x3 




- 1 ' 


-2 3(er)'-4-"4-"" + 3(e.-")V'' 


(44b) 


(44c) 


3 (4“^> 

The exact (cnmbersome) expression for the ratio can be obtained by plngging in Eqs. 


(70a), (70b) and (70c) into Eq. (44c). 


6.A. Universal aspects of the limiting distribution 


We now consider the c ^ oo properties of the limiting distribntion. Consider the ratio 
E[r^]/E[r]. To that end, we isolate the highest order term of c in the denominator and 
nnmerator and obtain 


E[r2] _ 2 2^^YBs + 0{cY 
E[r] “ 3 2^^cABs + 0{(A)' 


(45) 


We arrived at this expression by manipnlating the expressions for E[r] and E[r^] given by 
Eq. (60) and Eq. (69) respectively, that involved the terms in Eqs. (70a)- (70c). Therefore, 


-= - 

c-s-oo E[r] 3 


(46) 
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This limiting ratio is universal in the sense that it does not depend on fo and coincides with 
the answer obtained by integrating the DMPK distribution [5,17,27 


We will now compute the hrst two moments of the DMPK distribution in Eq. ([^. Let 
us suppose that the eigenvalues of the transfer matrix are log-uniformly distributed so that 


h{X)X = /tXjg i/g] for some small positive e such that e -C 1. Then Eq. (37) gives us 

= 1 


E [r] = Ak 


(A + iy 


;dX + 0{e) = 4k + 0(e), 


whereas Eq. (43) gives us 


E[r^] = 16 K 


A 16 

(y + 0(e) = —k + 0(e), 


so that 


E|t] 4 3 ^ 


and for e -C 1, we get the universal limiting ratio predicted Eq. (46). Thus in the c —)■ oo limit 


the DMPK distribution exhibits the same universal ratio of the hrst and second moments as 
the limiting distribution we have derived using random matrix theoretic arguments. 


The discussion in Section [6^ suggests that whenever the medium is ‘sparse’ in the sense 
that kM fV|ay/M —)■ c, we can expect to get a distribution of the form posited by the DMPK 
theory irrespective of the material properties of the individual scatterers. The natural next 
step in this line of inquiry is to analyze the large c asymptotics of the transmission coefficient 


distribution via its implicit characterization in Eq. (36) to answer hner questions about the 


existence of a density at A = 1 (equivalently r = 1) for all c G (0, cxd). We leave these for 
future work. 


6.B. Multiple-point-scatterer-per-layer scenarios that lead to same limiting distribution 

We now consider multiple-point-scatterer-per-layer scenarios that lead to the same limit¬ 
ing distribution - this will suggest a sparsity condition for the existence of the perfect 
transmission-supporting universal limiting transmission coefficient distribution. Consider the 
setting similar to that in Fig. except with k randomly placed point scatterers per layer. 
Then if D 3> r is large, we expect the S*!! and S 22 matrix to be approximately rank k, by 
neglecting the scatterer-scatterer interaction related terms. Consequently, we can model the 
empirical eigenvalue distribution of ■ Ti as 

V / j=i 
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Retracing the steps after Eq. (26), we observe that we arrive at the same limiting distribution 


encoded in Eq. (36) except now with kN\ay/M —)■ c and 


AW 


i=i 


Now, suppose that we are in the setting where the rank of the S*!! and S 22 matrices 
depends on M. Let us make this dependence explicit by denoting it as kM- Suppose that 


kM/M —>■ 0. Then, following the argument following Eq. (47), we will arrive at the same 


limiting distribution encoded in Eq. (36) whenever kj^ N\ 3 y/M —)■ c and 


AW 

Km , 


i=i 

Our analysis thus suggests the sparsity condition —)■ 0 and kM N\ay/M —)■ c for the 

emergence of the perfect transmission-supporting universal transmission coefficient distribu¬ 
tion. 


7. Numerical simulations 

To validate the predicted transmission coefficient distribution, we adopt the numerical sim¬ 


ulation protocol described in 13 . Specihcally, we compute the scattering matrices in Eq. 


(|^ via a spectrally accurate, T-matrix inspired integral equation solver that characterizes 
helds scattered from each cylinder in terms of their traces expanded in series of azimuthal 


harmonics. As in 13 , interactions between cylinders are modeled using 2D periodic Green’s 


functions. The method constitutes a generalization of that in 16 , in that it does not force 


cylinders in a unit cell to reside on a line but allows them to be freely distributed throughout 


the cell. As in 13 , all periodic Green’s functions/lattice sums are rapidly evaluated using 


a recursive Shank’s transform using the methods described in 22,23 . Our method exhibits 


exponential convergence in the number of azimuthal harmonics used in the description of 


the held scattered by each cylinder. As in 13 , in the numerical experiments below, care was 


taken to ensure 11-digit accuracy in the entries of the computed scattering matrices. 

We now describe how the simulations were performed. We generated a random scattering 
system with r = 0.05A,f' = 12.63A, D = 25.75A, and M = 51. The locations of the scatterers 
were selected randomly as described in Section]^ For a given A'^iay, the number of layers in 
the scattering system, we numerically compute the scattering matrices. We then compute 
the empirical transmission coefficient distribution over 200 Monte-Garlo trials and compare 
it to the analytically predicted transmission coefficient distribution obtained as a hxed point 


of Eq. (36) for c = N\^y/M and an appropriate choice of fo. 
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We first consider the setting where all the randomly placed cylinders have the same re¬ 
fractive index. Plugging in /^(t) = 5{t — 6) into Eq. (29) yields the expression 




+ 


z-d z-9 


-1 


z-1 


(48) 


For n = 2.08, we get a = 0.33 and 9 = 0.5. Plugging m. 9 = 0.5 into Eq. (48) and solving 


Eq. (36) yields the transmission coefficient as a function of c. Fig. shows the agreement 
between the physically rigorous empirical distribution and the analytically predicted dis¬ 
tribution. Note in particular, the agreement for c = 2 where the distribution is far from the 
characteristically bimodal DMPK distribution. 

We now consider the setting where with probability pi a cylinder has refractive index rii 
and with probability p 2 h has a refractive index n 2 . Plugging in fo{t) = pi5{t—9i)+p25{t — 92) 


into Eq. (29) yields the expression 




Pi 


z — 9i 


+ 


Pi 


01 


-1 


+ 


P2 


Z — 9r 


+ 


P2 


02 ^ 


z-1 


(49) 


For ni = 1.28 and n 2 = 2.89 we get cki = 0.05, 9i = 0.9, and 02 = 0.82, 6*2 = 0.1. 


Plugging in these values into Eq. (49) with pi = 0.8 and p 2 = 0.2 and solving Eq. (36) 


yields the transmission coefficient as a function of c. Fig. shows the agreement between 
the numerically obtained empirical distribution and the analytically predicted distribution. 
Note in particular, the agreement for c = 2 where the predicted distribution is supported on 
two intervals and agrees with the empirical results. 


Finally, we consider the setting corresponding to fd(t) = 


~^ 9 i<t< 92 - We generated 


02 — 01 

the scattering system by mapping each random realization of 0 to a random realization of 


the refractive index. Plugging this choice into Eq. (29) yields the expression 




02 — 01 


log 


01 


09 


+ 


02 


-01 1 , 
- + ^ log 


02-2 - 1 
9iZ — 1 


(50) 


The choice of 0i = 0.1 and 02 = 0.9 corresponds to a refractive index of rii = 2.89 (with 
tti = 0.82) and a refractive index n 2 = 1.28 (with 0^2 = 0.9). Plugging these values of 0i and 


02 into Eq. (50) and solving Eq. (36) yields the transmission coefficient as a function of c. 


Fig-i shows the agreement between the physically rigorous empirical distribution and the 
analytically predicted distribution. Note in particular, the agreement for c = 2 where the 
distribution is far from the characteristically bimodal DMPK distribution. 

Appendix [^contains some movies that shows the evolution of the transmission coefficient 
distribution with c for each of the three scenarios discussed. As expected, for large enough 
c the distribution eventually becomes characteristically bimodal as predicted by the DMPK 
theory. The behavior for small values of c is accurately predicted by our theory. 


17 

























For the three settings described above, we analytically compute E[r] from the associated 
fe via Eq. (40). The computation involves the normalizing factor B 2 , which for the three 
settings is given by 


B. 


nonrandom 


D atomic ^ 

B2 = Pi 


B. 


uniform 


= 1 - 


i-o 

1-^1 
1 + 01 
4 


+ P2 


log 


1-02 
1 + 02 
02 + 1 


+ 


02-01 ^V^l + V (01 + 1)(02 + 1) 


(51a) 

(51b) 

(51c) 


The closed-from expression for E[r^] is lengthy and therefore omitted here. It can be obtained 
using the calculations in Appendix Fig. § compares the empirical moments with the 
predicted moments and shows the good agreement for a range of values of c. 

Finally, we numerically validate the analytical prediction in Eq. (46). To that end, we 
generated a random scattering system with D = 197A,r = 0.11A,L = 3.4 x = 

430, 000, = 1.3, and M = 395. The locations of the scatterers were selected randomly and 

produced a system with I = 6.69A, where I is the average distance to the nearest scatterer. 
Let L denote the thickness of the scattering system we are interested in analyzing. We vary 
L from A to L and for each value of L we compute the scattering matrices associated with 
only the scatterers contained in the (0, L) portion of the (0, L) system we have generated. 
This construction ensures that the average density per “layer” of the medium is about the 
same. We computed the hrst and second moment of the empirical transmission coefficient 
distribution by averaging over 1700 random realizations of the scattering system and com¬ 
puted ratio as a function of c = M/L. Fig. shows that the empirical result validate our 
theoretical prediction. 
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Appendices 

A. Derivation of Eq. ( |10[ ) 

Here, we uncover the relationship between the singular value squared of S' 21 , r, and the 
singular value squared of T, A. Our derivation follows the approach in Section l.C.l]. 
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(a) C=2. (b) C=41. 



(d) C=41: Zoom-in. 


Fig. 3. The transmission coefficient distribntion for the setting where fe{t) = 
S{t — 6). The red line is the theoretical prediction - the histograms are from 
the physically rigorous simulation averaged over 100 trials. Note the agree¬ 
ment with theory in the C = 2 where the distribution is far from the DMPK 
distribution. 
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(a) C=2. 
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(b) C=33. 
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(c) C=2: Zoom-in. The predicted distribution is supported on two intervals. 



Fig. 4. The transmission coefficient distribntion for the setting where fe{t) = 
0.86{t — 0.9) + 0.26{t — 0.1). The red line is the theoretical prediction - the 
histograms are from the physically rigorous simulation averaged over 100 trials. 
Note the agreement with theory in the C = 2 where the distribution is far from 
the DMPK distribution. 
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(a) C=2. 


(b) C=25. 



X 


(c) C=2: Zoom-in. 



(d) C=25: Zoom-in. 


Fig. 5. The transmission coefficient distribntion for the setting where fe{t) = 

- - 7 r^ 0 i<t< 02 - The red line is the theoretical prediction - the histograms 

are from the physically rigorons simnlation averaged over 100 trials. Note the 
agreement with theory in the C = 2 where the distribntion is far from the 
DMPK distribntion. 
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Fig. 6. The first moment versus c for the settings corresponding to Fig. m Fig. 
I^and Fig. Irrespectively. The results of the physical simulations were averaged 
over 100 trials. 


A.A. Decomposition ofT^ ■ T 

Recall that the eigenvalues of ■ T equal the square singular values of T. Using Eq. @ 
and the fact that ■ S = I, we can express ■ T as 


rjnH rjn _ 


CH 

^21 

cH 
~ *^11 

. q-h . 

*^12 

cH 
' *^22 



cH 
■ *^22 


/ + 

25"- 

q-h . 
*^12 

51-2^ ■ 


—9 9“" ■ 9-1 
^^12 ^12 

■5n 


hi ^12 

'i-H 


-H 


S{2 


■Su -25"-R 


1- 

c c c 

*^21 “ *^22 ■ *^12 ■ *^11 

522 ■ 5r2' ' 


- 5 : 2 ' ■ 5n 

1_ 


25: 


11 ^12 ^12 
*^12 I 


12 


(52) 


In order to factorize the matrix on the right-hand side of Eq. (52) further, we first factorize 
the submatrices of the scattering matrix as 


521 = f/ ■ S ■ (53a) 

5n = F ■ U* ■ V/-S2 . U" (53b) 

5i 2 = F ■ 5ji ■ F = F ■ U* • S ■ (F ■ U*)^ (53c) 

522 = f/ ■ F ■ V/ - E2 . (F ■ U*f, (53d) 


where F = diag({e-i'^"}„) and £ [0,27r], and F represents the phase ambiguity between 
the singular spaces. Note that the factorizations in Eq. (53) satisfies power conservation, 
reciprocity and time-reversal symmetry. Substituting these into Eq. (52) yields the factor- 
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Ratio 
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Fig. 7. The ratio of E[r^]/E[r] as a function of c. The 2/3 line corresponds to 
the prediction in Eq. (46) for the large c limit of this ratio. 
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ization 


V ■ {2E-^ - I) ■ -2V ■ VI - ■ {F ■ V*)^ 

-2F ■ V* ■ V/- S2 . S-2 .1/^ F-V* ■ (2S-2 - /) . (F ■ V*)^ 


V 0 


2 S- 2 -/ -2VI - E 2 . S -2 


V 0 

0 F-V* 


-2VI - E 2 . S -2 2 S- 2 -/ 


0 F-V* 




/ 



=:E 


(55) 


Note that this factorization reveals that the eigenvalues of ■ T are exactly equal to the 
eigenvalues of S. 


A.B. Eigenvalues of a special block matrix 


Note that the matrix S on the right-hand side of Eq. (54) is of the form 


Di D 2 

D3 D4 


where Di = diag{{diViii), -^2 = diag({(i 2 ,j}i^i), -D 3 = diag({d 3 ,i},^J and D 4 = 
dia.g{{diV^i)- The eigenvalues z of this block matrix are the solutions of the character¬ 
istic equation 


det 


Di — zl D 2 

D 4 — zl 


0 . 


Since the eigenvalues will not be the same as di^i, Di — zl will be invertible; hence the 
characteristic equation can be rewritten as 


det{Di — zl) ■ det(Zi )4 — zl — ■ {Di — zl) ^ ■ D 2 ) = 0 . 


Equivalently, 

M M , d d \ ^ 

■ U (^ “ d^'"-z ) ^ n + di^idA,i - d2,id2„i] 

i=l i=l k 1.* / j=i 

Consequently, the 2M eigenvalues 2 ; are given by 

di i F dii ± V{dii + d^i)'^ — d{di idii — d2,id^V r • 1 nr 

z = -^- tor ? = 1 ,..., M. 


0 . 


(56) 


24 
























A.B.l. Relationship between the singular values of T and S '21 

Recall that Ti = af is the eigenvalue of and Aj is an eigenvalue of ■ T. We can apply 


Eq. (|56|) to the matrix S in Eq. (|55|) to obtain the eigenvalues of ■ T, which are given by 


^ (2/r. - 1) + (2/n - 1) ± y(4/T. - 2r - 4((2/r. - 1)" - 4(1 - t,)/^ 


= 2/ri-l±2,/l/T=-l/T, 


(58) 


Note that 1/Aj = 2/rj — 1=F — I/tj. This tells us that the singular values of the 

transfer matrix come in reciprocal pairs. Consequently, if we specify the singular values 
above one then the singular values below one are given by their reciprocal. 


B. Derivation of Eq. (40) for E[r 


From Eq. (39) the hrst moment is given as 

E[r| = - 4 a(-l). 


In order to evaluate this further, we are going to use what we have driven in Eq. (|33|), 

c) = exp ( -c • "I" 

V z{z + l)^ 

Using the relationship between ^ph{z) and ^h{z) in Eq. ( [^ , we get 

—1/ 


, 1 + ;^ 1 1 + ;^ ( e(l + ^)' 

i^) = —TTTT^ = — exp I C ■ 


Therefore, 


2 : i>h{z) 


E[r] = -4e;(-l) = -4. 


z{z + l) 


dzih i^)U=z 


(59) 


(60) 


where z* is a value that satishes ^/i(—1) = z* or ^^^(z*) = —1. We can easily check that 


= —0.5 by plugging it into Eq. (59). 


^(-0-5) = ^- 71 ^ exp ( c • I = -1 exp ( -4c • ^(-1)) . (61) 


-0.5 


-0.5(-0.5 + l) 


From Eq. (29), we can evaluate 'C(—1) as follows 

0.5 


e(-l) = - 


+ 


0.5 


-l-t -1-1/t - 1-1 

-0.5 -0.5f 


fe{t)dt 


t+1 t+1 


0.5 


fe{t)dt = 0 


(62) 

(63) 
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Therefore we confirm 


4‘(-0-5) = -lexp (-4c-?(-!)) = -1. 


(64) 


To complete Eq. (60), we need to evalnate dzi'^^{z)\z=-Q. 5 - Straightforward evalnation 
leads to the following result, 

2 


(^)|.2=-0.5 — —4 1 + C 


1 -t 
1 + t 




and we get 


E[Tl = -4a(-l) = -4^ 


S.C (x)|.=-o.5 1 ^ (PlC, 


C. Closed-form expression for ]E[r^] 
For notational convenience we define 


1 + t 


(65) 


Bn ■ = 
— 1 ' 


1 - t 
1 + t 
-1 


fe{t)dt 


■= Uh 


From Eq. (43) we have 


E[r^] = 16 




As we did for the first moment, we are going to express this in terms of the inverse 
of We are going to use the following elementary results from calculus 






-1'' 


3 fc"' 


C(-l) = 


-iv-i" 


- C‘ ik 




(66a) 

(66b) 

(66c) 

(66d) 

(67) 

function 

(68a) 

(68b) 


Substituting Eq. (|68a|) into Eq. (|67|) yields the expression 

T 


E[r^] = 16 


^C(-i) - ^e;:(-i) 


8 3(c^'? + 3e.- 

3 


C' a 




(69) 
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Note that and ^ can be expressed in terms of Bn, defined in Eq. (66a) as 


= -2\1 + cB 2 ) (70a) 

= -2\l + 2 cB2 + c^B^) (70b) 

^ = —2^ (3 + 6 ci ?2 + 6 (c + c^)B 4 + 2c?’B^ . (70c) 


The full (messy) expression for the second moment can be obtained by plugging Eqs. (70a), 
(70b) and ( 70c| ) into Eq. (69). 


D. Movies showing evolution of the transmission coefficient distribution with c 
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